package oneD.fem.algorithm;

/**
 * this is a replacement for BLAM in Windows for function which BLAM got error
 */

import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.jlinalg.IVectorReadOnly;

public class MyLinalg {

	/**
	 * perform y <- alpha*x + y
	 * 
	 * @param alpha
	 * @param x
	 * @param y
	 */
	public static void add(double alpha, IVectorReadOnly x, IVector y) {
		for (int i = 0; i < y.getSize(); i++)
			y.set(i, y.get(i) + alpha * x.get(i));
	}

	/**
	 * perform y <- alpha*x + y
	 * 
	 * @param alpha
	 * @param x
	 * @param y
	 */
	public static void add(double alpha, double[] x, double[] y) {
		for (int i = 0; i < y.length; i++)
			y[i] += alpha * x[i];
	}

	public static double dot(double[] x, double[] y) {
		double tmp = 0.0;
		for (int i = 0; i < x.length; i++)
			tmp += x[i] * y[i];
		return tmp;
	}

	public static void scale(double alpha, IVector x) {
		for (int i = 0; i < x.getSize(); i++)
			x.set(i, alpha * x.get(i));
	}

	public static void multiply(double alpha, int transA, IMatrix A, IVector x,
			double beta, IVector y) {
		if (transA == BLAM.NO_TRANSPOSE) {
			for (int i = 0; i < y.getSize(); i++) {
				double temp = 0;
				for (int j = 0; j < A.getColumnCount(); j++)
					temp += (A.get(i, j) * x.get(j));
				temp *= alpha;
				temp += (beta * y.get(i));
				y.set(i, temp);
			}
		}
	}

	public static void scale(double alpha, IMatrix M) {
		for (int i = 0; i < M.getRowCount(); i++)
			for (int j = 0; j < M.getColumnCount(); j++)
				M.set(i, j, alpha * M.get(i, j));
	}

	public static void scale(double alpha, double[] v) {
		for (int i = 0; i < v.length; i++)
			v[i] *= alpha;
	}

	public static void zero(double[] riGlobal) {
		for (int i = 0; i < riGlobal.length; i++)
			riGlobal[i] = 0;
	}

	public static double[] toArray(IVector v) {
		double[] a = new double[v.getSize()];
		for (int i = 0; i < v.getSize(); i++)
			a[i] = v.get(i);
		return a;
	}

	public static double norm2(double[] v) {
		double d = 0;
		for (int i = 0; i < v.length; i++)
			d += v[i] * v[i];
		return Math.sqrt(d);
	}

	public static void printVectorDouble(double[] b) {
		for (int i = 0; i < b.length; i++)
			System.out.print(b[i] + ",");
		System.out.println();
	}

	public static IMatrix inverseDiagonalMatrix(IMatrix M) {
		IMatrix Mi = (IMatrix) M.clone();
		for (int i = 0; i < Mi.getRowCount(); i++)
			Mi.set(i, i, 1 / Mi.get(i, i));
		return Mi;
	}

	public static void inverse2x2(IMatrix A) {
		double a = A.get(0, 0), b = A.get(0, 1), c = A.get(1, 0), d = A.get(1,
				1);
		double detA = (a * d - b * c);
		A.set(0, 0, d / detA);
		A.set(1, 1, a / detA);
		A.set(0, 1, -b / detA);
		A.set(1, 0, -c / detA);
	}

	public static void solve2x2(IMatrix A, IVector b) {
		double a1 = A.get(0, 0), b1 = A.get(0, 1);
		double a2 = A.get(1, 0), b2 = A.get(1, 1);
		double c1 = b.get(0), c2 = b.get(1);
		double detA = (a1 * b2 - a2 * b1);
		double x = (c1 * b2 - c2 * b1) / detA;
		double y = (a1 * c2 - a2 * c1) / detA;
		b.set(0, x);
		b.set(1, y);
	}

	/**
	 * solve the 2 unknowns equation a1*x+b1*y=c1; a2*x+b2*y=c2
	 * 
	 * @param a1
	 * @param b1
	 * @param c1
	 * @param a2
	 * @param b2
	 * @param c2
	 * @return double[2]{x,y}
	 */
	public static double[] solve2x2(double a1, double b1, double c1, double a2,
			double b2, double c2) {
		double detA = (a1 * b2 - a2 * b1);
		double x = (c1 * b2 - c2 * b1) / detA;
		double y = (a1 * c2 - a2 * c1) / detA;
		return new double[] { x, y };
	}

	/**
	 * compute determinent of 2x2 matrix
	 * 
	 * @param a1
	 * @param b1
	 * @param a2
	 * @param b2
	 * @return the determinant
	 */
	public static double det2x2(double a1, double b1, double a2, double b2) {
		return a1 * b2 - a2 * b1;
	}

	/**
	 * compute determinant of 3x3 matrix
	 * 
	 * @param a1
	 * @param b1
	 * @param c1
	 * @param a2
	 * @param b2
	 * @param c2
	 * @param a3
	 * @param b3
	 * @param c3
	 * @return the determinant
	 */
	public static double det3x3(double a1, double b1, double c1, double a2,
			double b2, double c2, double a3, double b3, double c3) {
		return a1 * (b2 * c3 - b3 * c2) - b1 * (a2 * c3 - a3 * c2) + c1
				* (a2 * b3 - a3 * b2);
	}

	/**
	 * copy x->y
	 * 
	 * @param x
	 * @param y
	 */
	public static void copy(double[] x, double[] y) {
		for (int i = 0; i < x.length; i++)
			y[i] = x[i];
	}
}
